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I. INTRODUCTION 



The study of the dynamics of suspensions of charged particles is interesting both because 
of the subtle physics underlying many electrokinetic phenomena and because of the practical 
relevance of such phenomena for the behavior of many synthetic and biological complex 
fluids In particular, electrokinetic effects can be used to control the transport of 

charged and uncharged molecules and colloids, using electrophoresis, electro-osmosis, and 
related phenomena j3j. As micro-fiuidic devices become ever more prevalent, there are an 
increasing number of applications of electro-viscous phenomena that can be exploited to 
selectively transport material in devices with mesoscopic dimensions j^. 

In virtually all cases of practical interest, electroviscous phenomena occur in confined 
systems of a rather complex geometry. This makes it virtually hopeless to apply purely 
analytical modeling techniques. But also from a molecular-simulation point of view electro- 
viscous effects present a formidable challenge. First of all, the systems under consideration 
always contain at least three components; namely a solvent plus two (oppositely charged) 
species. Then, there is the problem that the physical properties of the systems of interest 
are determined by a number of potentially different length scales (the ionic radius, the Bjer- 
rum length, the Debye-Hiickel screening length and the characteristic size of the channels in 
which transport takes place). As a result, fully atomistic modeling techniques become pro- 
hibitively expensive for all but the simplest problems. Conversely, standard discretization 
of the macroscopic transport equations is ill suited to deal with the statistical mechanics of 
charge distributions in ionic liquids, even apart from the fact that such techniques are often 
ill-equipped to deal with complex boundary conditions. 

In this context, application of mesoscopic ("coarse-grained") models to the study of elec- 
trokinetic phenomena in complex fluids seem to offer a powerful alternative approach. Such 
models can be formulated either by introducing effective forces with dissipative and random 
components, as in the case of dissipative particle dynamics (DPD) or by starting from 
simplified kinetic equations, as is the case with the lattice-Boltzmann method (LB). 

The problem with the DPD approach is that it necessarily introduces an additional length 
scale (the effective size of the charged particles). This size should be much smaller than 
the Debye screening length, because otherwise real charge-ordering effects are obscured by 
spurious structural correlations; hence, a proper separation of length scales may be difficult 
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to achieve. A Lattice-Boltzmann model for electroviscous effect was proposed by Warren p] . 
In this model, the densities of the (charged) solutes are treated as passive scalar fields. Forces 
on fluid element are mediated by these scalar fields. A different approach was followed in 
Ref. j?!, where solvent and solutes are treated on the same footing (namely as separate 
species). This method was then extended to couple the dynamics of charged colloids to 
that of the electrolyte solution. As we shall discuss below, both approaches have practical 
drawbacks that relate to the mixing of discrete and continuum descriptions. 

The LB model that we introduce below appears at first sight rather similar to the model 
proposed by Warren. However, the underlying philosophy is rather different. We propose to 
consider the fluxes between connected nodes as the basic physical quantities that determine 
the evolution of local densities. Such a formulation ensures local mass conservation, does 
not rely on fluxes or gradients computed at the lattice nodes (which constitutes a source of 
error in other models due to the need to approximate them on a lattice), and by choosing 
a symmetric formulation for the link fluxes in terms of the nodes that are affected, we 
can recover the proper equilibrium without spurious fluxes. Our model relies on an LB 
formulation for mixtures. Hence, the improvements of the formulation based on link fluxes 
will overcome some of the limitation of previous LB models for mixtures based on gradient 
expansions of a free energy jsj. 

The method described is very flexible, and in particular general boundary conditions are 
easily implemented. This feature makes also the proposed formulation attractive, since it 
avoids problems related to mass and charge conservation at fluid-solid interfaces, an artifact 
that has plagued previous LB implementations. It is then possible to model the dynamics 
of colloidal particles and polyelectrolytes in solution. The electrostatic interaction between 
them is derived from the charge distribution in the fluid. Hence, we do not need to assume 
any specific form for the interaction between charged colloids, or between monomers in 
a polyelectrolyte. Electro-osmosis, the sedimentation potential, electrophoresis or other 
electrokinetic phenomena can be easily treated within the model. In this paper we consider 
the first two to illustrate the capabilities of the method. 

The electrolyte is treated at the Poisson-Boltzmann level. We are not restricted to the 
linearized Debye-Hiickel regime and can study the electrokinetic effects at high charge den- 
sities, being only limited by ionic condensation (as happens for example in cylinders). The 
model we introduce will miss effects due to charge correlations. 
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The remainder of this paper is organized as follows. In Section 2 we describe the hy- 
drodynamics of fluid mixtures to set the general background. In Section 3 we describe the 
proposed numerical method and subsequently, in Section 4, we discuss how to model general 
solid interfaces within this lattice model. Section 5 focuses on the special case of interest to 
treat electrokinetic phenomena. In Section 6 we validate the method by analyzing different 
situations of interest, including electro-osmosis and sedimentation. 



II. HYDRODYNAMIC DESCRIPTION OF NON-IDEAL FLUID MIXTURES 

In some respects, the dynamics of electrolytes at hydrodynamic scales is analogous to that 
of multicomponent mixtures. The simplest electrolyte model consists of two ionic species and 
a neutral solvent. In order to provide the general framework for the description of electrolyte 
dynamics, we first briefly review the dynamics of mixtures on hydrodynamic length and time 
scales. As in all hydrodynamic descriptions, the starting point of any discussion are the laws 
of conservation of mass and momentum. 



A. Mass conservation 



Every species of the fluid mixture satisfies the usual mass conservation law: 

dpk 



dt 



+ V • pfcVfc = 0, (1) 



where is the velocity and pk the density distribution of the species labeled by k. The 
total density, p = J2k Pk, is also conserved, and satisfies an equation analogous to Eq. (dJ 
with respect to the barycentric velocity pv = J2kPk^k, which describes the evolution of a 
fluid element. If we refer the motion of all species to this common velocity, then Eq. (0) can 
be expressed as 

^^'+V-p,v = -V-j., (2) 



dt 

where we have introduced the relative current of species k, = Pfc(vfc— v), which accounts for 
all dynamical effects arising from the mismatch in velocities between the different species. 
On very short time scales, such currents are controlled by friction relaxation. However, 
for mixtures composed of molecular constituents (as is usually the case in electrolytes), the 
inertial time scale is extremely small; hence the relative current can be assumed proportional 
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to a thermodynamic driving force, which is proportional to the gradient of the chemical 
potential. As a result, the relative current of species i becomes diffusive and can be expressed 
as 3 

ik = -Yj Dikpk'^Pnk, (3) 

i 

where f5 is l/ksT, with ks the Boltzmann constant and 1/T the inverse temperature. (3nk = 
log pk + is the chemical potential decomposed in an ideal and excess part, while Dik 
corresponds to the diffusion coefficient that determines the flux of species i induced by spatial 
variations in the chemical potential of species k. For the sake of simplicity, we concentrate 
on the case where cross diffusion is neglected, and hence Di^ = Didik- By substituting 
the chemical potential in Eq. (jHl), we can then express mass conservation in the form of a 
set of convection-diffusion equations, expressing the two mechanisms that control density 
evolution for each species, 

dpk 



dt 



+ V ■ pkv = V ■ Dk [Vpfc + PfcV/J/i^] • (4) 



B. Momentum conservation 

Next, we consider momentum conservation. On the same length and time scales, momen- 
tum conservation implies that the barycentric velocity follows the Navier-Stokes equation, 
d 

— pv + V ■ pvv = r/VV + (V ■ v) - + F"^*, (5) 

tJL 

where rj and ^ are the shear and bulk fluid viscosities, respectively, while F^^* is the external 
force acting on a fluid element. The effect of the interactions among the different species 
enters as a net force expressed as the gradient of the local pressure, p. In the presence of 
spatial gradients, the pressure has in general a tensorial character, and can be derived from 
the free energy of the system. However, for ideal electrolytes, the local pressure can always 
be expressed as a scalar. Hence, for the sake of simplicity we will consider this situation in 
what follows. As a result, we only need to input the free energy of the mixture to determine 
both the pressure and chemical potentials. Speciflcally, if we know the free energy per unit 
volume /?/(r) = Ek Pk[log{A' Pk) - 1] + PT, then 

Pi^k = = \og pk + 13 ^f". 

= E pkP^'k -Pf = Y.(pk + pkPT) - pr. (6) 



where the free energy, the chemical potential and the pressure are position dependent. The 
first term of the pressure corresponds to the ideal-gas contribution, fip^'^ = J2kPk while the 
other two contain all the information of the interactions among the fiuid species. If there is 
one majority neutral component, which only contributes to the ideal part of the pressure, 
then the excess component of the pressure can be identified as the osmotic pressure of the 
mixture. In general, the pressure gradient follows from Gibbs-Duhem, 

/9Vp = Yl PkP^f^k = E ( Vpfc + PkPVpT) ■ (7) 

k k 

and acts as a force. We will use this interpretation in the LB implementation discussed in 
the next section. 

Using the last expression for the pressure gradient, the Navier-Stokes equation reads 

^pv + V ■ pvv = r]V\ + ^VV ■ V - Vp''^ - J2 Pk^Pp""" + F'^*. (8) 

k 

III. NUMERICAL LATTICE METHOD 

We propose a model that combines a description of momentum dynamics based on lattice- 
Boltzmann, with a numerical description of the convection-diffusion equation. Quantities 
are defined on the nodes of a lattice, r, and time evolves in discrete time steps. The lattice 
is prescribed by specifying its connectivity. The connections of each node are determined by 
specifying the set of allowed velocities, Cj, where the sub- index i runs over all the allowed 
velocities. Then, each node r is connected to the nodes r + Ci. 



A. Diffusion model 

For convenience, let us rewrite the convection-diffusion equation, Eq. (jU), in the form 

d 

^Pfc + V ■ pfcv = -V ■ jfc, (9) 

where the diffusive flux is 

jfc = -A(VpA:+PfcV/?pn- (10) 

For the sake of clarity, we discuss separately the change in density of the species k due to 
diffusion and to advection. The total change in time of the density is simply the sum of the 
two contributions. 
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1. Diffusion 



Let us assume for the time being that the mixture diffuses in a fluid at rest. Eq. Q then 
becomes 

d „ . 

g^Pk = - V ■ Jfc. 

Integrating both sides of this equation over a volume Vq and using the Green's formula 
Ivo ^ ■ = /ao j ■ ^dA, we obtain: 



where n is the outward unity vector normal to the surface, Aq, enclosing the volume Vq. 

As we have pointed out previously, we will consider densities defined on nodes of a lattice 
and the time evolution evolves at constant time steps. In this case, we can identify the 
volume Vq with the volume associated to that node, and Aq is related to the connectivity 
of the lattice nodes. Then, Eq. states that the change of the total number of particles 
enclosed in the volume corresponding to node r equals the sum of the outward fluxes. Such 
fluxes can only take place by mass transport to the neighboring nodes that are connected 
to the central node, according to the structure of the predetermined lattice connectivity. 
Hence, 



where nfc(r) is the number of particles of species k at node r, while jkii^) accounts for 
the fraction of particles of species k going to node r + Cj. If we consider the velocity 
moving opposite to i, i.e. Ci> = — Cj, we have jfcj(r) = —jki'ijc + Cj) because these fluxes 
are always defined considering that the particles move away from the reference node. This 
unambiguously show that the fluxes are related to the links joining the connected nodes, 
rather than being quantities defined on the nodes. 

It is worth noting that in the previous balance equation the relevant quantity is the 
number of particles of species k at node r, ^^(r), rather than its number density, Pki^). 
If we take the volume of a cell as our unit of volume, then pfc(r) = ^^(r). However, in 
the presence of solid boundaries this distinction may become relevant. The prefactor Aq in 
Eq. (fT^ is related to the geometrical structure of the lattice. Rather than connecting it 
directly with the area of the Wigner-Seitz cell that can be associated to node r, we derive 
its magnitude by computing how density diffuses to the neighboring nodes. In Section IVI Al 




(11) 




(12) 
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we will compute explicitly this geometric prefactor for a particular lattice. In the following, 
when referring to link mobility, we will use the symbol dk = DkAo. 

Using link fluxes to compute the variation of the densities of the different species avoids 
approximating the divergence on a lattice, a source of lattice artifacts, and the related 
potential spurious fluxes that may appear. Moreover, the use of these link fluxes also 
imposes locally mass conservation to machine accuracy, avoiding the errors caused by the 
discretization of the spatial gradient operator. We must still provide a prescription to 
implement the diffusive fluxes. These are in principle given by Eq. (fTUj) and involve spatial 
gradients between two neighboring lattice nodes. In equilibrium, nl^ ~ exp[—Pfi'^^] and, 
as a consequence, Eq. (|1(J|) predicts that all diffusive fluxes vanish. However, the direct 
implementation of Eqn. (jlUj) on a lattice will suffer from discretization errors that will result 
in small but noticeable spurious fluxes. To eliminate this effect, it is convenient to write the 
expression for the flux on a link as 



because, in this expression, the gradient becomes identically zero when the density distri- 
bution corresponds to its equilibrium form. This also holds for the discretized form to be 
discussed below. Consistent with the idea that the flux can be expressed in terms of link 
mass fluxes, we propose a symmetrized implementation of jki involving magnitudes defined 
in the two connected nodes, r and r + Cj. In particular, we write the flux of species k along 
the link c, as 



where Aj = |cj| = |cj/| is the distance between the two neighboring nodes. This symmetrized 
formulation ensures that, to machine accuracy, jfcj(r) = — Jfej'(r + Cj), and mass is conserved 
for the model elementary dynamic processes. Note that, based on the mass conservation 
expression, Eq. (fT^ . the global mass change of node r is the sum of the link fluxes, jki- 
Mass evolution in the diffusive limit is described only on the basis of mass flux divergence, 
as we have described. In general, the procedure developed based on link fluxes provides a 
consistent framework to obtain other gradients if needed. 
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(13) 
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2. Advection. 



Local density can also be altered due to advection if there is a local velocity of the fluid. 
If, for the time being, we disregard diffusion, the advection mechanisms can be written in 
the form. 



where v is the barycentric fluid velocity. In principle, the change in the number of particles 
could be computed on the basis of the advection along each link, in a way similar to Eq. (fT^ . 
However, as we will describe in the next section, the model we will introduce provides the 
velocity at each node, rather than the link velocity. In order to avoid numerical artifacts 
and spurious diffusion due to the interpolation to get such a link velocity, we propose an 
alternative implementation of the advection process. We still consider that nfc(r) give us 
the number of particles in a volume element centered around node r. Since we know the 
velocity of that node, v(r), in one step the node will virtually displace to r + v(r) . As 
a result, the volume associated to node r will intersect some neighboring cells of the real 
lattice (see FiglH). We then distribute the amount of particles into the intersected volumes 
proportionally to the intersected region. In Fig. ^ we depict in shadow the volumes that 
correspond to the fraction of the density that is transported in the new cells. The advantage 
of this approach is that it greatly reduces the spurious diffusion that usually results during 
advection in lattice models. To be more precise, even with the present method, advection will 
cause some spurious diffusion (proportional to the flow velocity). However, in Section IVI Al 
we show that, in practice, this effect is negligible. 

B. Lattice Boltzmann Method. 

In order to simulate the hydrodynamic flow of the fluid, we make use of the lattice- 
Boltzmann approach. This technique has been used extensively to model hydrodynamic 
flows in complex geometries It is equivalent to solving a discretized version of the 

Boltzmann equation with a linearized collision operator. This method describes the dynam- 
ics of a fluid in terms of the densities of particles that "live" on the nodes of a cubic lattice 
and have discrete velocities {ci}, where i labels the links between a lattice point r and its 
neighbors. The values of the velocities are chosen such that, in one time step, a particle 




(14) 
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moves along a link from one lattice node to its neighbor. In the Lattice-Boltzmann model, 
the unit of length is equal to the lattice spacing and the unit of time is equal to the time 
step. In addition, the unit of mass (or, equivalently, energy) is fixed by the requirement 
that, in the continuum limit, the transport equations for the lattice model approach the 
Navier-Stokes equation. This imposes a relation between the temperature and the speed of 
sound (see below Eg. ITHj). The central dynamic quantity in the Lattice-Boltzmann approach 
is the one-particle distribution function, /j(r,t), which describes the probability of having a 
particle at site r at time t with velocity Cj. The hydrodynamic variables are obtained as mo- 
ments of this distribution function over the lattice velocities, Cj; e.g. density and momentum 
can be obtained as 

p(r, t) = Yl /,(r, t), j(r, t) = p(r, t)v(r, t) = ^ c,/,(r, t), (15) 

i i 

respectively. 

In the presence of external forces, F, the evolution equation can be expressed as 

/,(r + c,,t + 1) = /,(r,t) + A, [/,(r,t) - /f (r,t)] + V. (16) 

where C[^] is a linear collision operator acting on which tends to relax the distribution 
function toward its equilibrium limit. Hence, one needs to specify the equilibrium distribu- 
tion as well as the collision operator. The collision operator ensures mass and momentum 
conservation (i.e. J2i '^ij = J2i CiCij = 0). Its eigenvalues also determine the viscosity of the 
fluid. The equilibrium distribution appearing in Eq. (fTHI) is that of an ideal gas. It can be 
shown that the Navier-Stokes equations are recovered keeping a low velocity expansion of 
the Maxwellian Q|, i.e. 

1 . 1 



eq 



p + ^Ci-} + ^pvv : (ciCi - cll) 



(17) 



where : is the double inner product, and the coefficients a\ depend on the geometry of 
the lattice, and are chosen to ensure that the anisotropy of the lattice does not affect the 
hydrodynamic behavior of the model, as well as ensuring that all the distribution functions 
are non negative. Moreover, is the speed of sound and its value depends on the values of 
the coefficients a\ but it is always smaller than unity (in lattice units). Finally, the term ipi 
accounts for the external force. It satisfies J^ii^i = ^ an d y',- Cj ijji = F. For a more detailed 
description of how to model the external force, see e.g. jll, 12|. 
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It can be shown 11| that in the hydrodynamic hmit one recovers the Navier-Stokes 
equation 



d 

— pv + V • pvv = ?7VVv + CVV ■ V - c^Vp + F. 



(18) 



Since the third term on the rhs is the pressure gradient for an ideal gas, if we fix the 
temperature such that ksT = c^, we then recover Eq. (jSj) for an ideal mixture. For non- 
ideal mixtures, we will introduce the missing contribution to the pressure gradients as a 
local external force, F. 

Introducing the mixture non-ideality as a local effective force implies that the fluid reacts 
with the appropriate susceptibility to applied external fields, although in the absence of 
spatial gradients the equihbrium distribution corresponds to that of an ideal gas. Since we 
are not concerned with local structure, the model can be regarded as an effective kinetic 
model, similar in structure to a linearized Vlasov equation. Hence, this approach differs from 
previous proposals which try to derive the hydrodynamics of non-ideal mixtures from kinetic 
models of mixtures |l3 or from a modification of the equilibrium distribution to recover the 
equilibrium pressure [8 1 . 

For a particular choice of the shear viscosity, r] = 1/6 in lattice units Q|, the general 
dynamic rule Eq. (fTHj) simplifies to. 



/.(r,t+l) 



p(r, t) + \ci ■ (j(r, t) + F) + -^pvv : {dd - c^l] 



(19) 



For the sake of convenience, we implement the model with this simplified updating rule. 
However, it is straightforward to implement the more general form that allows us to impose 
other values of the viscosity. 

The peculiarities of the non-ideality of the mixture enters through the forcing term (F)in 
Eq. (|19p. This forcing term can be decomposed into an external field and an interaction 
contributions, F = F^^* + F**"'. This interaction force, as previously described in Eq. ((7j), 
has the form F"*"' = J2k Pk^Pf^T- Using the same approach that we have used to model the 
convection-diffusion equation, we can determine the force acting on each link, Fj. Moreover, 
for the particular case where the diffusion matrix is diagonal. 



k L 



Jki 



Ukjr + Cj) - nfc(r) 

A,, 



(20) 



The advantage of using the force exerted on the links is that again, we keep a symmetric 
dependence on the neighboring nodes, and moreover -Fi(r) = — Fj'(r + Cj). Yet, in the 
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Lattice-Boltzmann update rule, we need the force acting on the node. This force can be 
obtained averaging the hnk forces, 

Ffir) = Y.a'c,^F,ir) a = x, y, z (21) 

i 

Let us now introduce an alternative way of treating the same systems. There are situa- 
tions, as is the case in electrolytes, where one of the components of the mixture is dominant, 
and plays the role of the solvent. In this case, we can single out this component, ps, and treat 
it separately from the rest. In particular, since ps ^ pk, we can approximate the overall 
density by the solvent density (p ~ ), and the overall momentum by the solvent momen- 
tum {pv = J2kPkVk — PsVs)- If we then relate the moments of the distribution function /j 
to the solvent density, i.e. J2i fi = Ps and J^i^ifi = Ps^s instead of Eqs. (fT3j) . we impose a 
constant solvent density in the incompressible regime. Hence, the rest of the components 
will need to compensate their densities to avoid any net local density variation. Although 
this incompressibility constraint is not exact, it may be a convenient approximation. From 
the point of view of the link force, Eq. (j2(Jj) . it has the computational advantage that one 
gets Fi = J2kjki/F>k and it reduces to the link diffusive flux previously computed, Eq. (fT^ . 
In this case the Navier-Stokes equation becomes 

^p,v, + V ■ p.v.v, = r/V V, + eVV ■ V, - clVps - keT J2 [Vp^ + Vp^] + F^"*, (22) 

k 

and by taking fc^T = c^, we recover an appropriate behavior when p^ ^ p^. 

The advantage of this approach is that densities of different species are dealt with on 
different footing, which may prove advantageous in certain applications, specially when 
dealing with boundary conditions that act differently on the solvent and solute, as it is the 
case if dealing with semipermeable membranes. Numerically, in this case there is a net force 
only when the density distribution deviates from its local equilibrium value, in contrast with 
the original method, where the density coming from the advection contribution balances the 
local force. This ensures an additional way to avoid spurious artifacts from the underlying 
lattice. 

IV. BOUNDARY CONDITIONS 

If the fluid mixture is confined between walls, or if colloids are added to the mixture, we 
need to specify how the densities and distribution function will interact with solid interfaces. 
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To account fully for such an interaction, we need to describe in turn how the distribution 
function behaves, how the particle number evolves, and how we estimate the interacting 
force at the surface. 

At a solid surface we expect hydrodynamic "stick" boundary conditions to apply. One 
way to impose these is to apply the so-called "bounce-back rule" on the links. However, the 
standard version of this procedure (see for example Ladd jn|) allows the fluid to leak into 
the solid. Although this leakage is usually innocuous, there are cases (a typical example 
being when electrostatics is part of the excess chemical potential) where this leakage may 
change the density of the solvent inside the solid, leading to a corresponding error in the 
pressure gradient . There exist alternative bounce-back rules that do not allow for any fluid 
leakage jl5l |. 

The formulation of our model in terms of link fluxes simplifies the implementation of 
boundary conditions for the fluxes of the different species densities, pk- Since the convection- 
diffusion equation involves only mass conservation, it is enough to impose that there is no 
net flux on any link that joins a fluid node and a solid node. We accomplish this by imposing 
that the diffusive flux jki = on such a link, and that the flux due to advection also vanishes. 
This second requirement is achieved by a kind of partial bounce-back move: the number of 
particles that would have been assigned to a solid node after advection is reflected back to 
its node of origin. 

The updating rule, both for the number densities of the convection-diffusion equations 
and for the lattice Boltzmann distribution function, requires the evaluation of gradients 
of chemical potentials. To this end, we need to specify the values of the excess chemical 
potentials on neighboring nodes, and those may involve the values of the fluid densities in 
contact with the solid wall. We consider that the relevant value of the density is that in 
contact with the wall, which is somewhere in between the fluid and the solid node. Such 
value can be obtained by requiring that it is consistent with the no-flux condition for the 
link flux of that species. The no flux condition is satisfied requiring (see Eq. |T!?|) ). 

nuiv + c) = nfe(r)e^K^('-+^')-'^^^W], (23) 

which should be understood as the extrapolation of the fluid density to ensure the absence 
of flux diffusion, and in general it is an implicit equation to obtain an estimate of the 
extrapolated number of particles, ?T,fc(r -|- Cj). Note that this fictitious extrapolated density 
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is a property of the link, not of the node. 

As we have mentioned in Section ITTTl the formulation based on the fluxes is based on the 
evolution of the number of particles contained in a given volume element. For the fluid nodes 
in the absence of solid interfaces the particle number is proportional to the number density. 
This is no longer the case close to a solid wall. This difference is pertinent because the 
excess chemical potential and the pressure are functions of the number density, pk- While 
for a wall at rest, one can still consider that the wall is equidistant from the nodes and Uk 
and Pk coincide, for a moving solid surface, the position of the solid boundary will change as 
it moves. In this case, a coefficient a that establishes how close the solid boundary is to the 
fluid node should be introduced. In the limiting case that the solid boundary is reaching the 
neighboring fluid node, the corresponding cell has a volume that is approximately half the 
volume of a usual cell, hence a = 1/2; in the opposite case when the solid surface reaches the 
solid node one gets accordingly a = 3/2. This coefficient then allows us to relate = apk- 
Although there exist different ways in which this coefficient may be computed, any smooth 
function that accounts for the volume change will be enough to avoid abrupt changes in the 
density when a fluid nodes is absorbed or created by the moving boundary. 

V. ELECTROKINETIC EQUATIONS. 

In the previous sections we have developed a model to simulate general non-ideal fluid 
mixtures. We will now analyze the special case in which the fluid mixture is an electrolyte. 
The simplest electrolyte model corresponds to a three-species mixture, two of them being 
the ionic species, p+ and p_ with charges z+e and Z-e, and the third one being the neutral 
solvent ps- e is the elementary charge, and 2+ and z_ are the valencies of the ions. The 
local charge can then be expressed as g(r) = e[zj^_p^{r) + z_p_{r)]. The simplest free-energy 
model corresponds to an ideal mixture in the absence of any local electric field, where we 
can write, 

f3f{r) = Y: Pk[^og{Alpk) - 1] + lf3q$ (24) 

k=±,s 

with $ being the electrostatic potential and the factor 1/2 avoids double counting. The 
chemical potential is then 

ppk = '^ogpk + pZk^, k = +,~, (3p, = \ogps 
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(25) 



The hydrodynamic evolution equations for this free energy model become, 



d r ^1 

■^Pk + V • pfcv = DkV ■ [Vpk + ezkPkV(3<l>\ (26) 
d 

— pv + V • pvv = r/VVv + e VV ■ pv - Vp + /? J] ezkPk^^- (27) 

We still need an additional equation that prescribes how the electrostatic potential is 
related to the local charge density. Since transport processes associated to mass and mo- 
mentum transfer in fluid mixtures are much slower than the propagation of electromagnetic 
waves, the electric field is completely determined by the Poisson equation 



B 



^kpk + Ps 

k=± 



(28) 



which has been expressed in terms of a dimensionless potential, $ = 6/5$, while Ib = 
(3e^ / (47re) is the Bjerrum length (the distance at which the electrostatic and the thermal 
energies are equal), with e the dielectric constant of the fluid. In the previous equation, ps 
stands for the charge density of the solid surfaces, if there are confining walls or moving 
suspended particles in the electrolyte. Obviously, a will be non-zero only on those solid 
surfaces. The Equations ()24|). (1^^). and are commonly referred to as the Electrokinetic 
equations. 

The electrostatic potential $ can be computed using standard techniques. Specifically, 
we have implemented a successive over-relaxation scheme (SOR)Q| as described in more 
detail in Ref . (Tf . The advantage of this model is that it does not presume a specific type of 
boundary condition, and can be easily generalized to deal with media of different dielectric 
constants. Although not as fast as other methods for solving the Poisson equation, it is 
adequate for our purposes because, once the local equilibrium charge profiles are achieved, 
the calculation of the disturbed electrostatic potential due to external forces is much less time 
consuming than the iteration part related to lattice-Boltzmann and convection-diffusion. 



VI. VALIDATION TESTS 



In order to validate the model that we introduced in the previous section, we compare 
its predictions against known results. In particular, we verify that the equilibrium charge 
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distribution is properly recovered on the lattice, and that out of equilibrium the different 
coupling mechanisms between fluid flow and charge inhomogeneities are properly accounted 
for. 

A. Effective diffusion 

As was pointed out below Eq. ()12|1 . the diffusion coefficient characterizing the discrete 
version of the diffusion equation is not the same as the link diffusion coefficient, dk but is 
related to it through a simple geometrical factor Aq that depends on the type lattice used. 
Aq can be evaluated as follows. Consider a situation where the transport of species k is 
purely diffusive. A density perturbation po, initially localized at node fq, will spread in 
one time step to the connected neighboring nodes. If the process is purely diffusive, we 
know the amplitude of the second moment of the density variation during this time step 
and X^i '^?p(ro + Ci,to + 1) = QDkPo = GAndkPri in a three dimensional cubic lattice. Let 
us consider for concreteness the D3Q18 lattice (isf, which is the lattice we used in our LB 
simulation. Since the link fluxes ji = d/^po/Ai, after one time step the density in each 
of the six nearest neighbors is dkpo, while the density in each of the other 12 connected 
nodes is dkPo/\/2. As a result, Y.i ^'iPi'^o + Ci, to + 1) = dk{6 + 12v^)po, which implies that 
Aq = 1 + 2\pl (or = dk{l + 2^(2)). Depending on the value of d^ it might happen 
that the total density transferred to the neighbors is larger than the initial density. For 
D3Q18 this gives us an upper bound for the input diffusion coefficient that ensures absolute 
stability, dk < 1/(6(1 + 2v^)) = 0.044. In practice, we find that for all cases that we have 
analyzed, numerical instabilities related to diffusion become relevant for values of the input 
diffusion coefficient dk > 0.05. In order to perform simulations at higher diffusivities, we 
need to modify the numerical scheme to simulate the diffusion equation. This instability 
can be overcome by introducing a multiple-timestep technique. To this end, we introduce a 
smaller diffusion coefficient dit = dk/Nn and iterate times the discrete diffusion equation, 
Eq. (HU, to advance the densities one time step. 

When applying this multiple timestep method to solve the lattice diffusion equation, one 
must compute carefully the force that should be applied to the distribution function /j at 
the end of the time step. In fact, F'^°' should be computed at all the intermediate steps. All 
these contributions should then be added to obtain the total force at the end of the iteration. 
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With this technique can vary the diffusion coefficient over several orders of magnitude. For 
example, in our simulations we could vary from Dk = 10"'^ to Dk = 6 (all in lattice 
units) . 

On top of the lattice effects on diffusion itself, advection can also induce spurious dif- 
fusion, because the lattice velocities do not coincide in general with the local velocity. As 
a consequence, a concentrated set of particles will spread over the lattice nodes, even if 
subject to a pure translational motion. Hence, only when the velocity is commensurate with 
the lattice spacing, both in direction and magnitude, will spurious diffusion be exactly zero. 
We must then quantify the amount of spurious diffusion. To this end, we consider an ideal 
binary mixture composed of a solvent with initial uniform density, Ps, and a solute with 
initial density pt- The mixture is contained between two parallel walls that are permeable 
to the solvent but impermeable to the solute. The fluid is moving with a uniform velocity 
V perpendicular to the walls. As a result of the impermeability of the walls to the solute, a 
steady state is reached, determined by the solvent density profile, pt{x), which satisfies 



Pt[x) = Poexp 



(29) 



where v is the fluid velocity, D* the effective diffusion coefficient, po the solvent distribution 
at contact with the wall located at xq- 

In Fig.Elwe show the effective diffusion coefficient measured by using Eq. ()29|1 as function 
of the fluid velocity for a range of values of the diffusion coefficient. We plot D*/Dq (where 
Dq is the diffusion coefficient for a quiescent fluid). In order to show that there exists an 
intrinsic advection-induced spurious diffusion, we plot in the inset of the same figure the 
difference between the effective and the input diffusion coefficient for many values of the 
input diffusion coefficient as function of the fluid velocity. Because all curves collapses, this 
graph shows that the diffusion coefficient induced by the advection depends only on the fluid 
velocity. We observe that the dependence on the (absolute value of) flow velocity is linear 
with slope 1/2. Following the procedure that we used above to compute the factor Aq, we 
can derive an expression for the advection-induced diffusion coefficient. In one dimension, 
a fraction vAt of the density p{x) is displaced to the next node, while a fraction (1 — v)At 
remains at the original node. The center of mass of the density is displaced by a factor 
vAt. Simple algebra then shows that the second moment of the density variation during a 
time step is < Aj >= v{l — v). The flow-induced diffusion coefficient in one dimension is 
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therefore D* = l/2v — l/2v'^. In three dimensions this expression is readily generahzed to 
yield 



By choosing a sufficiently low value of the flow velocity, and a sufficiently large value of Dq, 
we can largely suppress the effect of this advective diffusion. 

B. Electrolyte in a slit. 

Next, we consider a fluid confined between two parallel solid walls at rest, with a constant 
surface charge. The sht has a width L and the surface density charge is fixed to p(— L/2) = 



The space between the two slits is occupied by a solvent and counter-ions. In order to 
achieve global neutrality, the density of counter-ion is initially set to be uniformly distributed 
p(x) = -a/L, X e {-L/2,L/2}. 

The actual position of the hydrodynamic and electrostatic solid boundary cannot be re- 
solved within a lattice spacing. In the neutral case, for the viscosity and geometry considered 
the wall can be assumed to be halfway between two consecutive lattice nodes, as dictated 
by the bounce-back ruleQj. We will use this position as a reasonable approximation. In 
fact, the results we describe for a planar slit indicate that for a planar wall the electrostatic 
position of the wall can be taken as being midway between the boundary nodes. For a 
non-planar interface a separate calibration will be required. 

1. Equilibrium distribution of the counter-ion density 

In equilibrium, a uniform charge density on a flat wall will induce an inhomogeneous 
equilibrium density profile of the counter-ions. For this simple geometry, the charge-density 
profile of the counter ions is known analytically known (at least, at the Poisson-Boltzmann 
level) y, Q] for an arbitrary surface charge density: 




(30) 



p(L/2) = a/2. 



p{x) 



Po 



(31) 




where po = K'^/2ttIb, K is the solution of the transcendental equation 




(32) 
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which involves the wall charge density. Since we have an exact solution for the full Poisson- 
Boltzmann equation for arbitrary values of the wall charge, this geometry is a good case 
to analyze the limitations of the model dealing with large charges, i.e. beyond the linear 
Poisson-Boltzmann limit. For low surface charge densities, the linear regime is recovered by 
linearizing Eq. |221 and the parameter K becomes KunL = ^/A^dBO. 

In the opposite limit of high surface-charge density, K saturates at KgatL = vr. We can 
then quantify the deviation of the fluid from the linearized regime where the electrostatic 
interactions are small by analyzing the departure of KL from KunL. 

In Fig. OlA we show the equilibrium counter-ions distributions in both limits. In our 
simulations we fixed the Bjerrum length to be 0.4, the channel width to 20 lattice nodes, 
and we have varied the surface- charge density. In the plot we show the profiles for K/ Kun = 
1.01, 1.13 and 2.01, which correspond to a = 0.003125, 0.03125 and 0.3125 in dimensionless 
lattice units, respectively. The highest value of K is not far from the saturation value. The 
figure shows that, with the present method, we can indeed reproduce the correct counter-ion 
distribution, both in the linear and in the non-linear regime. In Fig. IHlB we compare the 
density profiles close to the wall in the non-linear regime for two different slit widths. The 
larger the surface charge the more localized the charge profile will be. The figure shows that 
increasing the resolution of the lattice does result in a small but significant improvement in 
the calculation of the charge distribution. Of course, the discrepancy would be greater for 
a more localized charge profile. In practice, only the computer resources (memory) will set 
an upper limit for the surface charge density that can be modeled reliably with the present 
scheme. 

2. Electro- osmotic flow 

Having verified that the model correctly reproduces the equilibrium behavior, we next 
turn to the calculation of flow caused by an external electric field. We apply a constant 
external electric field that is parallel to the slit, i?". This field causes hydro dynamic flow 
as it exerts a force on those fluid elements that carry a net charge. If we take y as the 
component along the walls and refer to x as the coordinate perpendicular to the walls, then, 
at the Poisson-Boltzmann level, the exact solution for the fluid flow in the steady state can 
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be written as 




log 



cos(^)J' 



cos{Kx) 



(33) 



where rj is the shear viscosity of the fluid. In our simulations, we model the constant electric 
field by taking into account the potential difference that it causes between neighboring lattice 
nodes (i.e. A^extiv) = E^Ay). 

Fig. m shows the computed electro-osmotic flow profile in a slit confined by hard walls with 
a charge density a = 0.003125 (in units of the elementary charge per square lattice unit). 
In the same figure, we also show the analytical solution (Eqn. (jHSl), with K/Kun = 1.01) 
that is exact in the Debye-Hiickel limit. Again, there is good agreement between theory and 
simulation. This suggests that the effect of electrostatic forces on the hydrodynamic flow is 
correctly taken into account in the simulations. 

C. Sedimentation velocity 

In the previous sections we have seen that the appropriate equilibrium charge distribution 
is reproduced both in the linear and non-linear regimes of the Poisson-Boltzmann equation, 
and that also a charge distribution induces the correct fluid profiles. We must still show 
that the opposite coupling works correctly, i.e. we must compute the hydrodynamic drag 
on a charged object in the absence of external electrical fields. 

To this end, we compute the sedimentation velocity of an array of charged spheres im- 
mersed in an electrolyte solution. In this case, the velocity of the colloidal particle induces 
a fluid flow that determines the steady charge distribution around the sphere. This charge 
distribution in turn affects the sedimentation velocity of the particle. Hence, all the different 
couplings between charge, electrostatic potential and fluid flow are present. Such a scenario 
las been analyzed previously with a different model [3] and analytically at infinite dilution 
isl ]. As a consequence, we can again check our simulations against known results. 

The system that we consider consists of a charged sphere of radius a in a three-dimensional 
box of size L. Because of periodic boundary conditions, this corresponds to a periodic 
array of spheres with volume fraction ip = {Ana^ / L^). In the simulation, we first allow the 
electrolyte to equilibrate with the particle at rest in the absence of external forces; hence 
the system develops its equilibrium double layer. Then, we apply the gravity as an external 
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body force applied to the fluid, i.e. we move in the system of reference of the colloid. In 
this way we avoid the problem of updating the particle's position due to its motion jioj ]. 
By forcing the colloid to be at rest, we will not conserve momentum, but by computing the 
mean fluid velocity in the steady state (which is reached on a timescale of order L'^p/rj), we 
can obtain the sedimentation velocity. 

We have flxed the Bjerrum length to Ib = 0.4 and the radius of the sphere to a = 4.5 in 
lattice units. We performed calculations for two different values of the solvent fluid density, 
Ps = 1, and ps = 20, while the density of the added salt pk was varied between 1.8 * 10^^ 
and 4 * 10"^ As we vary the salt concentration, we also change the Debye length from 3.3 
to 21. Since Ps ^ Pkj we have performed most calculations using the second version of our 
simulation scheme, as described at the end of Section 1111 Bl However, we also performed 
some simulations using the original model (taking the solvent density as the overall density). 
The only difference that we observe between the two implementations is a small variation 
in the numerical value of the sedimentation velocity. However, this difference already shows 
up for sedimentation of a neutral sphere. It is due to a small change in the fluid viscosity 
that is caused by a small difference in the overall fluid density in the two implementations. 
The valency of the macro- ion was chosen to be Z=10 which corresponds to the small charge 
limit. Although our computational scheme should also work outside the Debye-Hiickel limit, 
we restrict ourselves to this regime, because it is only in this limit that we can compare with 
existing analytical results. Speciflcally, Booth predicted that the sedimentation velocity, 
Uq{Z), of a weakly charged sphere of valency Z in the dilute limit can be expressed as jl^ 

^ = 1 - c,Z\ (34) 

where Uq is the sedimentation velocity of a neutral sphere, and C2 is a constant that can be 
computed analytically in the Debye Hiickel limit. For the simplifled situation of monovalent 
CO- and counter-ions, 2;+ = —z_ = 1, which have the same diffusivity, = D_ = D, the 
expression for C2 simplifles to 



where f{K,a) is a linear combination of exponential integral functions [3] and is a function 



of the inverse Debye length, k = A^^ = y4:TTlBJ2k ^kPk- We have checked that the sedi- 
mentation velocity scales as predicted with the viscosity. We have also verifled that we are 
indeed in the linear regime where the sedimentation velocity is proportional to the applied 
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gravitational field. In particular, for the two values of the density considered, ps-, the linear 
regime was obtained for forces per unit of volume such that the flow velocity never exceeded 
0.1 in lattice units. 

Fig.^lshows the sedimentation velocity of a weakly charged sphere {Z = 10) as a function 
of the inverse Debye screening length. As can be seen from the figure, the sedimentation 
velocities scales with the ionic diffusivity in the way predicted by Eq. ESI The inset in the 
same figure shows that this scaling breaks down at higher colloidal charges ( Z = 100), i.e. 
outside the range of validity of the linearized Poisson-Boltzmann description. 

Fig. ini shows the reduced sedimentation velocity {U^p{Z) /U^{Z = 0)) as a function of 
K,a for a range of volume fractions. As the volume fraction decreases, the curves approach 
Booth's infinite-dilution result, while the minimum sedimentation velocity moves toward 
the minimum value predicted by theory. In order to compare quantitatively the simula- 
tion results with Booth's theory, Eq. we must extrapolate the computed values for 
Uip{Z)/Uip{Z = 0) from the finite <y9-values of the simulations to the infinite-dilution limit, 
Uo{Z)/Uq{Z = 0). For neutral spheres Hashimoto has shown that that the sedimentation 
velocity converges very slowly to its infinite-dilution value, namely as j^] 

Z oj = 1 - 1-7601^^/^ + ^ + 0{^'), (36) 



Ladd has numerically verified this dependence |2J|. For charged spheres, due to the elec- 
trostatic screening, we still expect that the dominant (p dependence comes from excluded 
volume; previous results indicate that this is indeed the case 3|. When performing the dilute 
limit expansion, we therefore decided to single out the major volume fraction dependence 
by normalizing the simulation results with the Stokes drag coefficient, i.e. computing the 
low density limit of U^{Z)/Uo{0). As a result, it is reasonable to obtain the same functional 
dependence on (f as Hashimoto with a slightly different amplitude. Specifically, we expect 

^ = 1 - (1.7601 + 6)^^3 + 0(^2/3)^ (37) 

where e is much less than one. Eventually, the dilute limit is obtained by extrapolating 
Eq.EZIto = 0. 

In Fig. El we show the extrapolated sedimentation velocities for a particular value of Ka. 
The estimated error in the limiting sedimentation velocity is rather large. It could have been 
reduced by computing more values of the sedimentation velocity at low volume fractions. In 
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addition, there is some uncertainty in the value of the effective sphere radius. In the hght 
of these uncertainties, the agreement with the Booth hmit in Fig. IHl is gratifying. 



D. Absence of spurious fluxes 



We pointed out in Section UTTl that one of the incentives for developing the present model 
was to eliminate any mixing of continuous-space gradients and discretized gradient opera- 
tors. The reason is that the inevitable approximations associated with the discretization of 
gradient operators usually lead to the appearance of spurious mass and momentum fluxes, 
even in equilibrium. Such spurious fluxes are present in particular whenever there exist spa- 
tial inhomogeneities related for example to the presence of liquid interfaces. In the present 
approach, we only use lattice-gradient operators that have been constructed such that, in 
equilibrium, no flow can result. To demonstrate the effect that this has, we compare the 
present method with an existing "mixed" method. In particular, we consider a spherical 
colloid of radius a = 4.5, at rest in an electrolyte in a cubic box of diameter L = 20. The 
valency of the sphere is Z = 10 and the system as a whole is electrically neutral. In Fig. [7| 
we show the projection of the momentum flux in the equatorial plane of the sphere and 
compare these residual fluxes both for the model introduced in this paper and the model of 



Ref. 



Fig.[7|a shows that spurious currents, although small, are certainly not negligible in 
this case. Moreover, their magnitude is clearly correlated with the distance to the colloidal 
particle: the largest currents appear in the region where the spatial gradients are largest. 
For highly charged spheres (i.e. outside the linear Debye-Hiickel regime) these spurious 
fluxes will become larger. In contrast, in Fig. [7|b (present model), the spurious fluxes are at 
the level of machine precision. In fact, to make them visible at all, we had to multiply the 
momentum fluxes by a factor 10^^ relative to the old model. In other words: the residual 
fluxes are controlled by machine accuracy. Even at this level one cannot detect a correlation 
between the fluxes and the position of the sphere. We can conclude that the proposed model 
eliminates the appearance of spurious equilibrium fluxes. 
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VII. CONCLUSIONS AND DISCUSSION. 



We have introduced a new model to simulate the collective dynamics of non-ideal fluid 
mixtures, with special emphasis on its use to study electrokinetic phenomena. The method 
relies on a lattice- Boltzmann model, where the interactions are introduced as effective forces. 
In this respect, our model resembles a Vlasov kinetic model, as opposed to previous kinetic 
lattice models. In our approach the fluxes between neighboring lattice nodes are the fun- 
damental dynamical objects that couple external fields to both electrical conduction and 
hydrodynamic flow. 

As a result of the symmetric formulation of the flux between neighboring nodes we can 
impose strict local mass conservation. As a consequence, the present model is free of spurious 
boundary fluxes that plague all other lattice-Boltzmann models of fluid mixtures. Moreover, 
a link-based description has the additional advantage that boundary conditions are easily 
implemented. 

Secondly, by using a multi-step approach, we can vary ionic mobilities over many orders 
of magnitude. This feature of our model allows us to explore electroviscous effects over a 
wide range of Peclet numbers. We have shown that flow causes spurious advection-diffusion. 
However, this effect is well understood and can be made negligible in most practical cases. 

We have checked the performance of the model by studying equilibrium diffuse layers, 
showing that it is possible to recover both low and high charge density regimes. In the latter, 
the only limitation is related to computational resources, because a finer grid is required to 
resolve the narrower charge profiles that develop nearly highly charged walls. To test the 
coupling of electrostatics and fluid flow, we have computed the sedimentation velocity of 
a charged sphere. These simulations indicate that the existing theoretical predictions are 
reproduced in the low-charge, low-density limit. As the charge of the colloid is increased, 
the simulation results start to deviate from the theoretical predictions that apply in the 
linearized Poisson- Boltzmann regime. 

Even though in the present paper we have focused on electrostatic interactions and, in 
particular, we have not discussed molecular interactions that favor demixing, such interac- 
tions could also be incorporated in the present model. 
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FIG. 1: Density redistribution due to advection. To advect the charge of a given node (in this 
case, node number 5) in one time unit, we shift the whole cell with the local velocity vector of that 
node, {vxiVy). Next, we displace a fraction of density equal to the area of the cell that is now in 
the corresponding site. In the graph a fraction of the density equal to the shadowed rectangle area 
{vxVy) goes from cell 5 to cell 3, a fraction (1 — Vx)vy goes to cell 2, (1 — Vy)vx goes to cell 6, and 
(1 — Vx){l — Vy) stays at node 5. For the sake of clarity, the figure shows a two-dimensional flow. 
In practice, the analogous procedure is carried out in 3D. 
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FIG. 2: In the Lattice-Boltzmann model, advection causes some spurious diffusion. Tlie figure 
shows the computed effective diffusion, D* / Dq as function of the fluid velocity for the steady state 
described in Section Fill A 11 The curves are drawn for different diffusion coefficients at zero velocity: 
Dq = 0.38 (circles), Dq = 0.57 (squares), Dq = 0.76 (diamonds), and Dq = 1.34 (triangles). In the 
inset we show that the amount of diffusion induced by the flow does not depend on the equilibrium 
coefficient and has, for small velocities, a linear velocity dependence. Symbols are simulation result 
and the dashed line corresponds to the theoretical expression described in the text. 
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FIG. 3: Equilibrium distribution of the charge density of counter- ions (no added salt) in slit between 
two charged walls at a distance L. The abscissa measures the distance from the wall in units of 
L. The local density is expressed in units of the average charge density in the bulk: po = c/L. 
A) charge distributions for three values of the dimensionless parameter KL (see text): KL = 
0.553 (circles), KL = 1.57 (squares) and KL = 2.77 (diamonds). In the same figure, we have 
indicated the corresponding analytical results (Eq. I31|) (dashed curves) for a slit of width L = 20 
lattice spacings. Circles and squares correspond to the linear regime {K/Kun = 1.01 and 1.13 
respectively), while diamonds are close to the saturation limit (K/Knn = 2.01). B) The accuracy 
of the numerical solution for the charge profile can be improved by increasing the spatial resolution 
of the lattice, in this case from L=20 (diamonds) to L=40 (circles). Again, the analytical result is 
shown as a dashed curve. The curves in B correspond to the result for a highly charged surface, 
KL = 2.77 {K/Kun = 2m). 
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FIG. 4: Electro-osmotic flow profile in a slit of width L = 20 lattice spacings. The surface charge 
density, a = 0.003125 {K/Knn = 1-01), corresponds to the linear regime. The fiuid in between the 
slit contains only counter ions. The electric field is along the ?/-direction. It has a strength of 0.1 in 
units ksT/ (Ale), where Al is the lattice spacing and e is the elementary charge. The simulation 
results are compared to the theoretical prediction, Ea.f IHS)). shown as a dashed curve. 
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FIG. 5: Reduced sedimentation velocity of a periodic array of colloids of valency Z = 10 in 
an electrolyte as a function of Ka. The figure shows the results for two different values of the 
ionic diffusion coefficients. The curve for Dq^^ = 0.95 (circles) has been rescaled to the curve for 
Do = 0.19 (x) according to Eq. i.e. U{D) = [/(Z^J^^) * (dJ^V^o)- The superposition of the 
two curves shows that the scaling in obeyed. In the inset we also show the results for a colloid of 
valency Z = 100. However, in this high-charge regime the sedimentation velocity does not scale 
with the diffusion coefficient in the way predicted by the linearized theory. 
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FIG. 6: Sedimentation velocity of a periodic array of spheres of valency Z = W and hydrodynamic 
radius a = 4.3. The Bjerrum length Ib = 0.4 (in lattice units). The diffusion coefficient of both 
positive and negative ions is set to D = 0.19. We compare simulation results for finite volume 
fractions, namely 0.0416 (squares), 0.0123 (diamonds), 0.00521 (triangles), and 0.00267 (circles) 
against the Booth theory, which is valid at infinite dilution (dashed curve). For ku = 0.5 we also 
show the estimated value of the sedimentation velocity at infinite dilution (see text). The point 
corresponds to the extrapolation of the law Eq. (|37|) . Within the estimated error, the extrapolated 
simulation results agree with the predictions of ref. [l^ . 
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FIG. 7: Illustration of suppression of spurious boundary currents in present LB model. In the 
figure we compare the apparent currents in equilibrium for two models: figure a) gives the results 
for the model described in ref. 9| , figure b) shows the results for the present model. In both cases 
we consider a colloidal sphere of radius 2.5 in a system with a diameter L = 20 lattice spacings. As 
there are no external forces acting on the system and the colloid is not moving, the fiuid is supposed 
to be at rest. The figure shows the measured projection of the momentum fiux in the equatorial 
plane of the colloid. In figure a), spurious currents are apparent close to the particle surface. The 
spurious currents in case b) are much smaller than in case a). In fact, to make them visible at all, 
they have been scaled up by a factor 10^^ with respect to case a). This is an indication that the 
spurious currents in case b) have been suppressed down to machine accuracy. 
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